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Waiting-time statistics are generated from the Olami-Feder-Christensen model and shown to 
mimic some aspects of real seismicity. Preliminary analysis of the model data implies a recently 
proposed universal scaling law for the distribution in seismicity may be due to a mixing between 
aftershocks and uncorrelated event pairs, thus having limited application. Earthquake catalog data 
is also presented to support the argument. 
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Earthquake systems display complicated statistics in 
both time and space, the specific mechanisms for which 
remain unclear. This is despite some apparently sim- 
ple power-laws describing them. The simplest, best 
known, and probably most significant of these laws are 
the Gutenberg-Richter (GR) [l| law for the probabil- 
ity of energy release, P(E) oc E~ b with b ~ 1, and 
Omori's Law for aftershocks Q which states that the 
rate of seismic occurrence decreases as R oc (t + c)~ p af- 
ter a large event, with p » 1. Many approaches have 
been considered to account for such statistics, includ- 
ing dynamical spring-block systems, cellular automata, 
and more generally the statistical frameworks of Self- 
Organised-Criticality (SOC) and the Epidemic Type Af- 
tershock Sequence model (ETAS). In this paper we are 
concerned with a dynamical model that is generally ac- 
cepted to be almost, if not precisely, SOC. The Olami- 
Feder-Christensen (OFC) model |3] is perhaps the most 
widely studied model of non-conservative self-organised 
criticality (if indeed it is technically SOC) (eg. 0, |]|), 
and has long been known to reproduce the GR law. More 
recently it has also been shown to qualitatively exhibit 
Omori's law as well as other features of seismic after- 
shocks dGJ . 

Correlations in time between earthquakes have long 
been recognised in the form of Omori clusters, but it is 
only recently that significant progress concerning corre- 
lation between the clusters has been made. The distribu- 
tion of waiting times for a set of events, which we'll write 
D(T), is often used to examine correlations in time. It 
may depend on a number of quantities, which will be 
discussed later. It is well known that completely uncor- 
related events, describing a Poisson process, result in a 
pure exponential decay for D(T). Variations from such 
indicate temporal correlations, and have long been looked 
for in earthquake catalogs. A short-time power law with 
exponent around -1 for D(T) is usually taken as related to 
Omori's law for the rate, while there is also evidence for 
longer term correlations amongst large events, eg. 0,0, 
in some cases, a faster decaying power law has been found 

Following , Corral has proposed a "universal" scal- 
ing law where D(T) is the same all over the world, dif- 



fering only according to the seismic rate, R, in the re- 
gion, so long as the rate is it is stationary. Specifically, 
the law states that D(T) = R x f(T/R) where / is a 
common function holds the world over. In [rl|. Corral 
calculated distributions from the NEIC catalog 0] for 
many regions and periods where the rate was approxi- 
mately stationary, and plotted them on a log-log graph 
where the axes had been rescaled by the rate. He then 
fit a single generalised gamma distribution to all the col- 
lapsed curves, specifically f(fi) = C i"L T e^^^ - 1 with 
7 = 0.67 ± 0.05, 5 = 1.05 ± 0.05, a =1.64 ± 0.15. This 
function is essentially a decreasing power-law of expo- 
nent 1 — 7 giving way to an exponential decay at longer 
times. This fit can be seen in Fig.l of each of 0, fl3j |. 
There is a deviation from / at short times, which is put 
down to rate fluctuations caused by aftershocks. It can 
be shown that a single Omori sequence should give rise 
to a Weibull distribution of waiting times . Although 
to the best of our knowledge this hasn't been confirmed 
in real events, we will use the result here. 

In a way, all the work concerning scaling laws started 
from a paper by Ito, ^5|, where the behaviour of a com- 
posite D(T) constructed from dividing southern Califor- 
nia into a grid of smaller regions was mimicked by the 
Bak-Sneppen model [lfjj|. In this way it is natural that 
we should look at a similar model of SOC to continue the 
work. We will show that the stationary rate data used by 
Corral may be well fit using the OFC model with certain 
parameters. Qualitative similarities are also found for 
more general behaviour. The model distributions show 
some complicated behaviour, but most can be described 
by three sections- a short time power law of exponent 
approaching —1, a transition period, and a rapid decay, 
fitting qualitatively with earthquake behaviour. We show 
evidence that pairs of events in the distribution may be 
separated roughly into those within a single aftershock 
sequence, and those that aren't correlated. We also find 
long time correlations in the model in the form of a steep 
power law, similar to that reported in 01- We usually 
find and exponent of around -2. 

The OFC model has been shown Q to be analagous 
to a 2-dimensional slider-block model (eg.^^) °f a 
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fault, and is described as a "continuous time cellular- 
automaton". A 2D lattice of side length L is defined 
representing an array of sliding blocks connected to each 
other by springs. The blocks are imagined to be driven 
uniformly by an overlying plate connecting to the top of 
the blocks by springs. The relative stiffness of the springs 
between each block and the vertical springs is represented 
by the only other parameter of the model, a, and governs 
the conservation of stress when a block slips. A value of 
a = 0.25 corresponds to the unphysical case of complete 
stress conservation, and a = corresponds to complete 
dissipation, that is the blocks are not connected to each 
other at all. 

The model works due to the inhomogcncity introduced 
by the boundaries. If not for this, as in the case of pe- 
riodic boundaries, the blocks synchronise into a state 
where all events involve only one block. When any 
amount of inhomogeneity is introduced, an ordered state 
extrudes from there, until the whole lattice is covered in 
larger patches, where each block in a patch has a simi- 
lar stress to it's neighbours. The patches are essential in 
understanding the model. They grow with distance from 
the boundaries, but do so much quicker with larger a, 
so that their average size is strongly dependent on both 
L and a. The state where the whole lattice is made up 
of such patches we'll call the critical state, although the 
technical correctness of the term is debated. It is in this 
state that the GR and Omori's laws are obeyed. Each 
patch relaxes in a succession of events, and the relax- 
ation of a single large patch appears to be the source of 
a foreshock/mainshock/aftershock sequence 0. 

For our work, data was only taken after each lattice 
had reached the critical state, as easily checked by dis- 
playing the lattice using colour to denote the stress on 
each block (as in the pictures in 0]). We have used dou- 
ble precision floating point values in all simulations for 
which data is displayed, although we have also considered 
single precision in light of recent studies > an d our re- 
sults are largely unmodified. We generated return-time 
data for lattices using open boundaries up to L=1024 
with a = 0.2, and smaller lattices using a as low as 0.05. 
We only present results for open boundaries here, as we 
found that free boundary conditions took much longer to 
simulate. This was due to the average event size being 
larger in the free case, which we put down to a longer 
range influence of the boundaries. A detailed analysis of 
mechanisms will be given in a another paper |2pj , but the 
behaviour of D(T) is qualitatively the same. Changing 
to free BCs and increasing L and M c will give a similar 
shape to using open boundaries. 

When taking data from the model we use the same 
scheme as Corral. We define a number of bins whose 
size increases exponentially toward longer times. When 
running a simulation, waiting times are placed in the ap- 
propriate bins, and the final count in each bin is divided 
by the bin width and the total number of events to get 
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FIG. 1: Results for a = 0.1, L = 256 and a = 0.2, L = 64 
lattices with various magnitude cutoffs. Fhe a — 0.2 curves 
are divided by 10 5 for clarity. Mc for these curves refers to the 
cutoff magnitude defined as the logarithm with base 2 of the 
number of slipped blocks. The inset data was calculated using 
Mc=8 and 11 for each lattice respectively. The line included 
is a power law with exponent -2 



the average probability density for that bin. To describe 
magnitude in the model, we use the logarithm of the 
number of blocks that slipped to the base 2. This is ap- 
propriate because limitations of size in the model result 
in a smaller magnitude range than seen for earthquakes. 

We can first make some general comments on the re- 
sults as plotted with logarithmic axes. The rate of oc- 
currence for any curve has the effect of moving a curve 
left or right. If events always occur in the same relative 
pattern, changing the unit of time should be equivalent 
to increasing the rate. Increasing L, a and decreasing 
M c all have the effect of increasing the rate. The precise 
value of M c is often unimportant because of the scale in- 
variance implied by the G-R law. The effect on the rate 
is clear: increasing L means a bigger area, and hence 
a higher rate. High a means less dissipation, so more 
stress remains in the lattice and more toppling occurs. 
Also, the general effects of L and a may be explained in 
a simple manner. The size and number of the patches 
is the important point here- increasing a increases the 
average patch size and hence the number of patches de- 
creases, as indicated by a smaller b value for the GR law. 
Decreasing L has a related effect, although the effect on 
the b value isn't as strong 

Fig. n shows various M c for two contrasting lattices. 
The L=64, a — 0.2 lattice consists mainly of a single 
large patch, whereas the L=256, a = 0.1 lattice has 
many. There are a few features to note in this graph. 
First, we look at the short time, hence low M c behaviour. 
There is a clear short-time power law in all three a = 0.1 
curves, with exponents of ~ —0.9. The a = 0.2 lattice 
displays an exponent w —0.5 when all events are consid- 
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FIG. 2: Results for the interior L/2 x L/2 section of lattices 
with a = 0.1, L = 256, and a = 0.2, L = 256. The a = 0.2 
curves are divided by 10 6 for clarity. Mc refers to the cutoff 
magnitude defined as the logarithm with base 2 of the number 
of slipped blocks. We are looking here for a decay in the 
distribution caused by aftershock sequences alone. 



ered. When larger events are looked at only, this part of 
the curve begins to flatten out. We find the exponents 
to be consistent with the exponent for Omori's law in 
these lattices, using the definition d=0 as described in 
0. This is consistent with the expected Weibull distri- 
bution of waiting times calculated from individual Omori 
sequences. This distribution should have a short-time 
power law corresponding to the Omori exponent |l4| . 
followed by a stretched exponential decay. If each after- 
shock sequence is completely defined by a single patch, 
then this is definitely the case. We will look for such a 
decay shortly. 

We can see that it is appropriate to divide the graph 
region into two areas, divided at about T w T a = 1 — 4a. 
Curves passing this point show a dramatic change in be- 
haviour. After this point, the curves in Fig^ become 
sharply oscillating decreasing functions, generally follow- 
ing a power law behaviour with exponent around -2. The 
inset shows this behaviour more clearly. The a = 0.2 
curve is much clearer. There are actually two periods of 
oscillation here. The first is « T a — 0.2 and the second 
is closer to 0.4 This corresponds to the average period 
in the middle of the lattice and at the boundaries. As 
discussed in j^, T a is the average period of each block 
in the interior in the limit of infinite lattice size. The 
blocks directly next to the boundaries have only three 
neighbours (for open BCs) and hence their average pe- 
riod should be 1 — 3a, in the middle of a side in the limit 
of a large lattice. As we do not use infinite lattices, we 
expect the period to be slightly less than this. While pe- 
riodic behaviour is quite clear in the a = 0.2 lattice, the 
a = 0.1 lattice shows dramatic drops at a long period in- 
terspersed with smaller scale oscillations. We think this 



may be due to a complicated beating behaviour due to 
many different periods. The oscillation occurs for large 
events because the patches change on a slow scale, and 
will relax in a similar manner on subsequent cycles. Gen- 
eral model behaviour such as this will be discussed in a 
future work. 

In Fig. ^the a = 0.1 £=256, and Mc=5 curve shows a 
very slight indent at T rj 10 -3 , although it is hard to see 
in that figure. The curves in Fig. 0where calculated only 
considering events initiated in the interior L/2 x L/2 part 
of the lattice. The idea is to capture the behaviour of the 
larger patches far from the boundary. In these we can see 
that the short time power law leads into a faster, non- 
exponential decay without interference from the periodic 
behavior at T a . Although it is reminiscent of a power 
law, we consider it more likely to be related to the sum 
and interference of separate Weibull distributions. 

Next we treat the model data in the same way as Cor- 
ral treated seismic data [n). We include data from 5 
different lattices. The first three are chosen to scale sim- 
ilarly to Corral's data. These lattices consist of only a 
single patch or so, while not being so small that they 
are dominated by boundary events. They show a low 
proportion of their events in the short-time power law. 
The other two are more typical, and consist of many 
patches. These show a large part power law behaviour. 
We haven't included curves which continue significantly 
past T a , as the behaviour produced prevents them from 
collapsing onto a single curve. Finally, we also include 
some curves from a lattice dominated by the boundarys. 
The a = 0.2, L — 32 lattice can't form patches prop- 
erly, and correlations dissapear. The exponential decay 
representing Poisson behaviour provides a reasonable fit, 
although perfect fit is impossible as the driving speed 
implies an inherent time-scale. 

As Corral's included data is chosen for it's stationary 
rate, it only offers a small insight into to thee variety of 
seismic behaviour. To get a fuller picture, we have cal- 
culated waiting time distributions in Fig^ for some re- 
gions of the world chosen for the goodness of their statis- 
tics, over periods for which they are complete. For this 
we used the combined cat alog of the Advanced National 
Seismic Network (ANSS), [23j. This catalog includes the 
NEIC catalog used by Corral. Duplicate events are re- 
moved according to the method described on their web- 
site. It should offer a more complete source than the 
NEIC catalog alone. 

The first curves are for events anywhere in the world. 
We can see that for lower cutoffs the scaling law is very 
well satisfied, but for larger cutoffs the short-time power 
law behaviour increases thus deviating at short times, 
while there is also a deviation from the exponential decay 
at longer times. The curves calculated from the Southern 
California and Alaskan regions of the ANSS catalog look 
quite similar to some of the curves in Fig. [21 The be- 
haviour of these catalogs we've found to be representative 
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FIG. 3: Various distributions plotted with axes rescaled by 
the rate. The top group are from lattices that give a good fit 
to the gamma function. These are L = 16 a = 0.1, L — 32 
a — 0.15, and L = 128 a = 0.2. Larger lattices or smaller 
a yields behaviour such as in the L=256 curves. Smaller 
L or larger a yields random behaviour. Subsequent sets are 
divided by multiples of 10 3 for clarity. Also shown for each set 
is f(n), or an exponential decay for the final set representing 
a Poissonian result. 




FIG. 4: Various distributions calculated from the ANSS 
worldwide catalog. The regions and periods shown are the 
world 1973-2004, Southern California 1984-2004, and Alaska 
1971-2004. The periods are chosen based on the catalog sub- 
mission details page of the ANSS website, (2^1 . The regions 
are also defined there. The Southern Californian and Alaskan 
curves are divided by 10 3 and 10 6 respectively. 



of many highly active regions of the world. 

To understand the scaling behaviour, we first note that 
the Poisson distribution is described solely by the rate, 
and gives pure exponential waiting-time distribution. Be- 
cause it depends only on the rate, scaling the axes it's 
plotted against in this way results simply in e~ x every 
time. In the model curves that scale well, the majority 
of events are in the exponential tail, so it is no surprise 



that they move together. As the number of aftershocks 
becomes significant, the power-law behaviour becomes 
significant resulting in poorer scaling. The model curves 
that scale well have the common property of consisting 
of only a few patches, and thus they have similar relative 
aftershock activity. The right amount of deviation from 
uncorrelated behaviour seems to be required to fit the 
gamma function, and it is not a general rule. 

As for earthquakes, we must work with a smaller num- 
ber of events, and definite conclusions are more difficult. 
We have shown however that when general behaviour is 
looked at, without isolating specific periods of station- 
arity, / is the exception rather than the rule as in the 
model. Further, when we look at the ANSS catalog data 
we see that the / is only well fit when a cutoff magnitude 
of about 5.5 is used. Lower gives a more random appear- 
ing curve, while higher gives more aftershock behaviour. 
This implies that the goodness of fit may be due to ei- 
ther incompleteness in the catalog, or simply a lack of 
correlations resulting in more randomness, when smaller 
events are considered. 

We have presented waiting-time statistics generated 
from the OFC model with open boundary conditions, 
showing that the general shape and scaling behaviour 
of real earthquake distributions is reproduced. It ap- 
pears that the generated distributions are formed from 
two regimes and so universal behaviour is not found. The 
reason for the scaling and the gamma function fit are then 
a majority of Poisson-like events, and a transition region 
between the regimes which approximates a slow decay 
power-law. Given the available evidence, it seems likely 
that real events scale for the same reason- when after- 
shock activity is low, that is the rate is approximately 
constant, most events follow a Poisson distribution. Also 
of note, we find that the waiting times of many large 
events describe an w — 2 exponent power law decay as 
has sometimes been found earthquakes. It is tempting to 
hypothesise that earthquakes may show this behaviour 
for a similar reason: that large areas of the earth remain 
in an organised state that changes only slowly with time. 
In this characteristic time governed by the slow 

driving would be apparent. In the model, it is the time 
taken for a block to go from zero stress to the threshold 
influenced only by the slow driving. For earthquakes, it 
could be the time taken for the friction in a single fault 
to build to it's slipping value, without being influenced 
by others. 

We should note that for our comparisons to real events, 
we are interpreting the model as representing a whole 
earthquake system, although it was originally intended to 
represent a single fault. This is however in line with the 
treatment in |2 1| on the fractal distribution of epicenters. 
In this case, the stress transfer may correspond roughly 
with the transfer of stress between sections of fault, al- 
though the analogy is somewhat lacking. The OFC model 
as it stands is a very unrealistic model of an entire seismic 
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system for a number of reasons, the main one being that 
real systems have no well defined boundaries, and are in- 
homogeneous across their area. The value of our results 
is perhaps not in revealing the fundamental behaviour 
of seismicity, but a rather more humble exploration of 
the mechanisms that may go in to constructing the often 
confusing waiting-time distribution. Further work will is 
required to determine if a real link to seismicity exists. 
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